NonParametric Kernel Density Estimation and Local Regressions with Mathematica

Introduction to the use of NonParametrix.m, a Mathematica Package.

Bernard Gress

University of California at Riverside

downloadable at: http://student.ucr.edu/~gressb01/index.html

NonParametrix.m, the Package

NonParametrix.m is a Mathematica Package that provides many of the basic (as well as a few advanced) functions often used in nonparametric econometrics and statistics, as described in, for example, Pagan and Ullah (1999), Silverman (1986), or Härdle (1989).

The functions included in this freely downloadable Package are:

NPKDense -

Provides nonparametric density estimation of multivariate data with user defined kernels.  Includes very fast routines for standard kernels such as the Normal, Epanechinikov, and the Uniform, for example.

NPRegress -

Fits an arbitrary order, nonparametric local polynomial regression estimation on multivariate data of any dimension, with user-defined or pre-defined kernels.  Also estimates slopes, regression residuals, confidence intervals, and can output an Interpolating function of the fits to use elsewhere.  Also includes local polynomial generalized least squares (GLS) routines.

CrossValidatedH -

Finds the optimal window-width ('h') for multivariate nonparametric regressions and density estimation.  Includes mean square error, root mean square error, median absolute error, mean absolute percentage error, and median prediction error loss functions.  Allows the user to select many of the parameters of the cross-validation routine.

Other Assorted Functions -

SilvermanH - Calculates the well-known asymptotic kernel window-width for data of arbitrary dimensions.  FastNPDensity - Provides very quick and dirty univariate density estimation when the user is mostly interested in quick visualization of data.  MultipleKDensityPlot - Provides very quick and dirty estimation of multiple densities, color coded for easy data visualization.

Needs["NonParametrix`"]

NPKDense

?NPKDense

Univariate Data

With no Options set then estimation is done at each data point, using the Silverman Rule-of-Thumb window-width:

ListPlot[Sort[Estimates /. NPKDense[unidata] ], 
        PlotJoined -> True];

NPKDense - Univariate Data (con't)

User defined Kernel

The user can supply his own Kernel where the x- y- z variables represent the (Subscript[ψ, i]); the variables must of course be mentioned in the domain declaration.  In this case the Kernel is the left-half of a Normal distribution:

f = (2/(Sqrt[2 π])) (E^(-((x^2)/2))); 
    domain[f] = {x, -∞, 0};

Plot[f, {x, -3, 0}, PlotRange -> All];

ListPlot[
        Sort[Estimates /. 
            NPKDense[unidata, Kernel -> f, 
              PointsToEstimate -> Chop[Range[- .2,  .2,  .002]]]], 
        PlotJoined -> True];

NPKDense

Bivariate Data

With user-supplied window-widths, estimation range, and number of points.  With MakeSquare turned on the output can be easily plot with ListPlot3D:

Timing[test = 
          Estimates /. 
            NPKDense[data2, H -> { .006,  .006}, 
              PointsToEstimate -> {- .1,  .1, 44}, 
              MakeSquare -> True];]

ListPlot3D[test, PlotRange -> All];

NPKDense - Bivariate Data (con't)

Bivariate Data

With a user-supplied Kernel, window-widths and number of points to plot, and the default estimation range:

f  =  (1/(x y σ (Sqrt[2 π]))) 
        Exp[-((((Log[x] - μ))^2)/(2 (σ^2)))] 
        Exp[-((((Log[y] - μ))^2)/(2 (σ^2)))]
domain[
          f]  =  {{x, 0, ∞} , {y, 
            0, ∞}};

Timing[test = 
          Estimates /. 
            NPKDense[data2, H -> { .014,  .017}, 
              Kernel -> f /. {μ -> 0, σ -> 1}, 
              PointsToEstimate -> 77, 
              MakeSquare -> True];]

ListPlot3D[test[[2]], PlotRange -> All];

ListContourPlot[Chop[test[[2]]]];

NPKDense

Trivariate Data

Another user-supplied, multivariate kernel, just to make things difficult:

g = 
      Simplify[PDF[
          MultinormalDistribution[{0, 0, 0}, #], {x, y, 
            z}]]
domain[
        g] = {{x, -11, ∞} , {y, -∞, 
          3}, {z, -5, ∞}}

Checking that we have an appropriate H matrix:

hh = {{ .01,  .001, 0}, {0.001,  .011, 0.002}, {0, 
          0.002,  .009}}
N[Eigenvalues[hh]]

Timing[test = 
          Estimates /. 
            NPKDense[data3, H -> hh, Kernel -> g, 
              PointsToEstimate -> {-1.3, 1.1, 15}, 
              MakeSquare -> True] ;]

ListContourPlot3D[test[[2]], Lighting -> False, 
        Axes -> True, 
        Contours -> {0.001, 0.01, 0.1, 1.0, 10. , 100. }, 
        ContourStyle -> {{RGBColor[0.4,  .4,  .8]}, {RGBColor[
                0.1,  .7,  .2]}, {RGBColor[ .8,  .7,  .1]}, {RGBColor[ .7,  .4, 0]}, {RGBColor[ .2,  .3,  .5]}}];

NPKDense - Trivariate Data

Ok, let's do a hard one with 40×40×40 = 64,000 points!:

Here we can see a contour plot,in 3D:The Joint Density of Stock Returns and its First Two Lags.  This is comparable to the Bivariate, 2-D contour plot of the distribution shown previously.

gg = 
        ListContourPlot3D[test[[2]], Lighting -> True, 
          Axes -> True, Contours -> {0.0001, 0.01, 1, 10}, 
          ContourStyle -> {{RGBColor[ .8,  .7,  .1]}, {RGBColor[
                  0.1,  .7,  .2]}, {RGBColor[0.4,  .4,  .8]}, {RGBColor[
                  1,  .4, 0]}, {RGBColor[ .2,  .3,  .5]}}];

click-and-drag on the images to rotate them in real-time... Quazar!

NPRegress -

?NPRegress

NPRegress -

The univariate version and arbitrary order case of Pagan and Ullah, pg 93 is:

\(Minimize\+{m, β}\ \(∑\+\(i = 1\)\%N\((\((y\_i - 
                m - \(β\_1\) \((x\_i - 
                      x)\) - \(β\_2\) \((x\_i - x)\)\^2 - … - \(β\_P\) \((x\_i - x)\)\^P)\) 
            K \((\(x\_i - x\)\/h)\))\)\)\)

which can be achieved by the weighted least-squares regression of

((((∑_(i = 1)^N) (Subscript[z, i]) (Subscript[K, i]) 
                ((Subscript[z, i])^′)))^(-1)) (∑_(i = 1)^N) 
              (Subscript[z, i]) (Subscript[K, i]) (Subscript[y, i])

where

(Subscript[z, i])=(({{1}, {((Subscript[x, i]) - x)}, {((((Subscript[x, i]) - x))^2)}, {⋮}, {((((Subscript[x, i]) - x))^P)}}))

and

(Subscript[K, i]) = K ((((Subscript[x, i]) - x)/h)).

Extending Pagan and Ullah, pg 93 to the multivariate and arbitrary order case:

\(Minimize\+{m, β}\ \(∑\+\(i = 1\)\%N\(\((\(\((\(y\_i - 
                          m - \(β\_\(1, 1\)\) \((x\_\(1, i\) - 
                                x\_\(1, j\))\) - \(β\_\(1, 
                                2\)\) \((x\_\(1, i\) - x\_\(1, j\))\)\^2 - … - \(β\_\(1, 
                                P\)\) \((x\_\(1, i\) - x\_\(1, j\))\)\^P -  ... \) - \(β\_\(k, 1\)\) \((x\_\(k, i\) - 
                            x\_\(k, j\))\) … - \(β\_\(k, 
                            P\)\) \((x\_\(k, i\) - x\_\(k, j\))\)\^P)\)\^2\) 
              K \((\(x\_\(1, i\) - x\_\(1, j\)\)\/h\_1, …, \(x\_\(k, i\) - x\_\(k, j\)\)\/h\_k)\))\)\(.\)\)\)\)

which can be achieved by the weighted least-squares regression of

((((Subscript[Z, i]) (Subscript[K, k, i]) 
              ((Subscript[Z, i])^′)))^(-1)) (Subscript[Z, i]) (Subscript[K, k, i]) 
      (Subscript[y, i])

where

(Underscript[((Subscript[Z, i])), ((P×k) + 1)×N])=(({{1, 1, ⋯, 1}, {((Subscript[x, 1, 1]) - (Subscript[x, 1, j])), ((Subscript[x, 2, 2]) - (Subscript[x, 2, j])), ⋯, ((Subscript[x, k, N]) - (Subscript[x, k, j]))}, {((((Subscript[x, 1, 1]) - (Subscript[x, 1, j])))^2), ((((Subscript[x, 2, 2]) - 
                      (Subscript[x, 2, j])))^2), ⋯, ((((Subscript[x, k, N]) - (Subscript[x, k, j])))^2)}, {⋮, ⋮, ⋮, ⋮}, {((((Subscript[x, 1, 1]) - (Subscript[x, 1, j])))^P), ((((Subscript[x, 2, 2]) - 
                      (Subscript[x, 2, j])))^P), ⋯, ((((Subscript[x, k, N]) - (Subscript[x, k, j])))^P)}, {((Subscript[x, 2, 1]) - (Subscript[x, 2, j])), ((Subscript[x, 2, 2]) - (Subscript[x, 2, j])), ⋯, ((Subscript[x, 2, N]) - (Subscript[x, 2, j]))}, {((((Subscript[x, 2, 1]) - (Subscript[x, 2, j])))^2), ((((Subscript[x, 2, 2]) - 
                      (Subscript[x, 2, j])))^2), ⋯, ((((Subscript[x, 2, N]) - (Subscript[x, 2, j])))^2)}, {⋮, ⋮, ⋮, ⋮}, {((((Subscript[x, 2, 1]) - (Subscript[x, 2, j])))^P), ((((Subscript[x, 2, 2]) - 
                      (Subscript[x, 2, j])))^P), ⋯, ((((Subscript[x, 2, N]) - (Subscript[x, 2, j])))^P)}, {⋮, ⋮, ⋮, ⋮}, {((Subscript[x, k, 1]) - (Subscript[x, k, j])), ((Subscript[x, k, 2]) - (Subscript[x, k, j])), ⋯, ((Subscript[x, k, N]) - (Subscript[x, k, j]))}, {((((Subscript[x, k, 1]) - (Subscript[x, k, j])))^2), ((((Subscript[x, k, 2]) - 
                      (Subscript[x, k, j])))^2), ⋯, ((((Subscript[x, k, N]) - (Subscript[x, k, j])))^2)}, {⋮, ⋮, ⋮, ⋮}, {((((Subscript[x, k, 1]) - (Subscript[x, k, j])))^P), ((((Subscript[x, k, 2]) - 
                      (Subscript[x, k, j])))^P), ⋯, ((((Subscript[x, k, N]) - (Subscript[x, k, j])))^P)}}))

and

(Underscript[((Subscript[K, k, i])), (N×N)])=(({{\(K \((\(x\_\(1, 1\) - x\_\(1, j\)\)\/h\_1, …, \(x\_\(k, 1\) - x\_\(k, j\)\)\/h\_k)\)\), 0, ⋯, 0}, {0, \(K \((\(x\_\(1, 2\) - x\_\(1, j\)\)\/h\_1, …, \(x\_\(k, 2\) - x\_\(k, j\)\)\/h\_k)\)\), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, \(K \((\(x\_\(1, N\) - x\_\(1, j\)\)\/h\_1, …, \(x\_\(k, N\) - x\_\(k, j\)\)\/h\_k)\)\)}})).

where (Subscript[x, k, j]) is the (k^th) dimension of the point at which local estimation is being performed.

NPRegress - Univariate Data

With no Options set the default is local-linear 
(Order = 1), with Silverman 'h':

MultipleListPlot[Last[Transpose[data1]], 
        Sort[Transpose[NPFits /. NPRegress[Last[Transpose[data1]]]]], 
        PlotJoined -> {False, True}];

NPRegress - Univariate Data (con't)

Order -> 0 ⟺ Nadaraya-Watson
and User window-width:

MultipleListPlot[Last[Transpose[data1]], 
        Sort[Transpose[
            NPFits /. NPRegress[Last[Transpose[data1]], H -> 4, 
                Order -> 0]]], 
        PlotJoined -> {False, True}];

NPRegress - Univariate Data (con't)

Epanechinikov kernel, with output of residuals and slopes:

results = 
      NPRegress[data1, H -> 0.6, Kernel -> Epakinikov, 
        Order -> 1, Residuals -> True, Slopes -> True]; 
    MultipleListPlot[data1, Transpose[NPFits /. results], 
      Transpose[NPResiduals /. results], 
      PlotJoined -> {False, True, False}];
ListPlot[
        Sort[Transpose[{First[Transpose[data1]], 
              First[Last[NPSlopes /. results]]}]], 
        PlotJoined -> True];

NPRegress - Univariate Data (con't)

Output of Interpolating Slopes Function:

results = 
        NPRegress[data1, H -> 0.6, Kernel -> Normal, 
          Slopes -> "Interpolated"];
MultipleListPlot[data1, Sort[Transpose[NPFits /. results]], 
        PlotJoined -> {False, True}];

testfunction = SlopeFunction /. results
Flatten[testfunction[First[Transpose[data1]]]]

NPRegress - Univariate Data (con't)

Output of the Interpolated Slopes for Order = 0.

In this case there are only Interpolated Slopes, and the Slopes -> "Interpolated" statement is redundant:

results = 
        NPRegress[data1, H -> 0.6, Kernel -> Epak, Order -> 0, 
          Slopes -> "Interpolated"];
MultipleListPlot[data1, Sort[Transpose[NPFits /. results]], 
        Sort[Transpose[{First[NPSlopes /. results], 
              First[Last[NPSlopes /. results]]}]], 
        PlotJoined -> {False, True, True}];

NPRegress - Univariate Data (con't)

Output of Confidence Intervals:

results = 
        NPRegress[data1, H -> 0.6, 
          ConfidenceIntervals -> 0.99];
MultipleListPlot[data1, Sort[Transpose[NPFits /. results]], 
        Sort[Transpose[{First[NPFits /. results], 
              First[NPConfidenceIntervals /. results]}]], 
        Sort[Transpose[{First[NPFits /. results], 
              Last[NPConfidenceIntervals /. results]}]], 
        PlotJoined -> {False, True, True, True}, 
        PlotStyle -> {Automatic, {Hue[0.1], 
              Thickness[ .009]}, {Hue[0.6], Thickness[ .006], 
              AbsoluteDashing[{Dot}]}, {Hue[0.6], Thickness[ .006], 
              AbsoluteDashing[{Dot}]}}, 
        SymbolShape -> {PlotSymbol[Diamond], None, None, 
            None}];

NPRegress - Univariate Data (con't)

Automatic Output of Confidence Intervals:

results = 
        NPRegress[data1, H -> 0.6, 
          ConfidenceIntervals -> {3, 0.958}, 
          OutputGraph -> True];
Show[NPGraph /. results, 
        DisplayFunction -> $DisplayFunction];

NPRegress - Bivariate Data

Order -> 0 ⟺ Nadaraya-Watson

Silverman 'h' automatically chosen...

Timing[results = 
          NPRegress[Flatten[data2, 1], 
            Order -> 0];]
ListPlot3D[
        Table[(NPFits /. results)[[2, dd (i - 1) + j]], {i, 
            dd}, {j, dd}]];

NPRegress - Bivariate Data (con't)

Confidence Intervals

The default values are used for the other parameters like H, Order, and PointsToEstimate.

results = 
        NPRegress[Flatten[data2, 1], H -> {0.9, 0.4}, Order -> 1, 
          ConfidenceIntervals -> 0.9999];
cc = 
        ListPlot3D[
          NonParametrix`Private`makeSquare [
            Flatten[Last[(NPFits /. results)]]], 
          Mesh -> False];

NPRegress - Bivariate Data (con't)

User-defined Kernel

The user can define any kernel he wishes to smooth with:

g = (E^(-((x1^2)/2))) (E^(-((y1^2)/2)));       domain[
        g] = {{x1, -3, 3}, {y1, -∞, 2}};

Timing[
      results = 
        NPFits /. NPRegress[Flatten[data2, 1], Kernel -> g2b, 
            H -> { .2,  .6}, 
            Order -> 1]; 
ListPlot3D[
        Table[results[[2, dd (i - 1) + j]], {i, dd}, {j, 
            dd}]];]

NPRegress - LLGLS

Local GLS Regressions

(Subscript[y, i]) = m((Subscript[x, i])) + (Subscript[ε, i])

    =m(
              x) + ((Subscript[x, i]) - x) β(x) + (Subscript[ε, i])

    =α(
            x) + (Subscript[x, i]) β(x) + (Subscript[ε, i])

where

α(x) = 
      m(x) - x β(x).

Therefore:

(Subscript[y, i]) = (Subscript[z, i]) δ(
            x) + (Subscript[ε, i])

where

(Subscript[z, i])=[1  (Subscript[x, i])]   and    δ(x)=(({{α(x)}, {β(x)}}))

or

y = 
      Z δ(x) + ε,

where Z is a matrix generated by:

(Subscript[z, i]) = [1,  (Subscript[x, i])].

To get the Local-Linear GLS estimator for δ, we minimize

\(TraditionalForm\`\(\((y - 
                Z\ \(δ(
                    x)\))\)\^′\) \(\(\(Ω\^*\)\^\(-1\)\)(
          y - Z\ \(δ(
                x)\))\) = \(ε\^′\) \(\(Ω\^*\)\^\(-1\)\) ε\)

\(TraditionalForm\`\(δ\&^\)(
        x) = \(\((\(Z\^′\) \(\(Ω\^*\)\^\(-1\)\) 
                Z)\)\^\(-1\)\) \(Z\^′\) \(\(Ω\^*\)\^\(-1\)\) \(\(y\)\(.\)\)\)

\(TraditionalForm\`\(Ω\^*\)\^\(-1\) = \(\@\(K(
              x)\)\) \(Ω\^\(-1\)\) \@\(K(x)\)\)

\(\(β\&^\_LLGLS\) \((x)\) = \(\((\(X\^′\) 
                K \(\((x)\)\^\(1\/2\)\) \(Ω\^\(-1\)\) 
                K \(\((x)\)\^\(1\/2\)\) X)\)\^\(-1\)\) \(X\^′\) 
        K \(\((x)\)\^\(1\/2\)\) \(Ω\^\(-1\)\) 
        K \(\((x)\)\^\(1\/2\)\) y\)

NPRegress - Benchmarks

Competing Applications: XploRe, N

({{Comparisons of Timings}, { (on Win2K, Athlon XP,  1.6  MHz)}, {All times in Seconds}, {({{ 
Package,  , 10  K Points

                    Univariate

                   Order  =  0, 10  K Points

                    Univariate

                   Order  =  1, 10  K Points

                    Univariate

                   Order  =  2,  }, {NPRegress -, \(\(\ \ \ \ \ \ \)\(Fast

                Local\ Poly

                \ \ Routine

                \)\), 
52.7, \(\(
\)\(60.2\)\), 
76.0,  }, { ,    Arbitrary

                Order/Kernel

                    Routine, 
319.7, \(\(
\)\(695.3\)\), 
1177.4,  }, { , 
Convolution

                    Routine, 
?, \(\(
\)\(NA\)\), 
NA,  }, {N,  , ?, ?, ?,  }, {XPloRe,  , ?, ?, ?,  }})}})

CrossValidatedH

Cross Validation searches for the window-widths that minimize some loss function in the prediction of some out-of-sample observations with a non-parametric regression, or, in the case of density data, the prediction of other datapoints

?CrossValidatedH

"CrossValidatedH[X, Y, Options] for the estimation of 'h' for non-parametric regressions, \!\(\*
StyleBox[\"or\",\nFontVariations->{\"Underline\"->True}]\) CrossValidatedH[Data, Options] for the estimation of non-parametric density estimations where X is a vector or a (K×N) matrix and Y is a Length-N vector, \!\(\*
StyleBox[\"or\",\nFontVariations->{\"Underline\"->True}]\) where Data is an (N×(K+1)) matrix.  Finds the optimal smoothing parameter \!\(\(h\^*\)\) based on the average out-of-sample predictive accuracy.  The search is over the range specified by 'HRange', considering the predictive accuracy over a 'Leaveout' number of dropped observations.  If 'Leaveout' is between 0 and 1, then it is interpreted as a percentage; otherwise it is the number of observations to drop.\\n\\nOptions are:\\nOrder -> \!\(\*
StyleBox[\"0\",\nFontVariations->{\"Underline\"->True}]\)|integer;\\nHRange -> \!\(\*
StyleBox[\"Automatic\",\nFontVariations->{\"Underline\"->True}]\)|(K×2) matrix of Mins and Maxs between which to estimate the grid of points;\\nHGridDimensions -> \!\(\*
StyleBox[\"10\",\nFontVariations->{\"Underline\"->True}]\)|integer specifying the resolution of the grid over which minimization of the LossFunction is performed;\\nLossFunction -> \"\!\(\*
StyleBox[\"MSE\",\nFontVariations->{\"Underline\"->True}]\)\"|\"MedianSE\"|\"MAE\"|\"RMSE\"|\"MAPE\"|\"MPE\" (mean squared error, median squared error, mean absolute error, root mean squared error, mean absolute percentage error, and median prediction error);\\nKernel -> \!\(\*
StyleBox[\"Normal\",\nFontVariations->{\"Underline\"->True}]\)|all kernels available to NPRegress;\\nLeaveOut -> \!\(\*
StyleBox[\"1\",\nFontVariations->{\"Underline\"->True}]\)|an integer number of observations to drop out-of-sample when measuring the LossFunction;\\nOutputLosses -> \!\(\*
StyleBox[\"True\",\nFontVariations->{\"Underline\"->True}]\)|False;\\nGraphLosses -> \!\(\*
StyleBox[\"False\",\nFontVariations->{\"Underline\"->True}]\)|True.  Use 'Show[CVGraph/.results, DisplayFunction->$DisplayFunction]' to display output;\\nReport -> False|True, outputs a little report summarizing the search for the optimal h;\\nMinMethod -> \"\!\(\*
StyleBox[\"Simple\",\nFontVariations->{\"Underline\"->True}]\)\"|\"Interpolation\"|\"NP\" - determines by what method the LossFunction will be minimized: by simple minima, by using a non-parametric interpolation;\\nSamples -> \!\(\*
StyleBox[\"15\",\nFontVariations->{\"Underline\"->True}]\)|All|integer - allows out-of-sample measurement of loss to be done over fewer samples than 'All' (i.e. N, the data size);\\nSearchPoints -> \!\(\*
StyleBox[\"Automatic\",\nFontVariations->{\"Underline\"->True}]\)|integer - Used in the NMinimize routine to determine the number of random beginning points.\\nInterpolationOrder -> \!\(\*
StyleBox[\"3\",\nFontVariations->{\"Underline\"->True}]\);\\n\\nOutputs are: {HStar, CVGraph, Losses, LossMinima}. 'Report' is included when it is set to 'True.  'HRange' is included whenever it is calculated automatically."

CrossValidatedH  - (con't)

Univariate Data

There are a wide variety of parameters that the user can fine-tune...

kernel = 
      Normal; {cvh, losses, g, hrange} = {HStar, Losses, CVGraph, HRange} /. 
        CrossValidatedH[data[[1]], 
          data[[2]], 
          LossFunction -> "MSE", Order -> 1, 
          Kernel -> kernel, LeaveOut -> 10, 
          GraphLosses -> True, OutputLosses -> True, 
          Samples -> 30, MinMethod -> "Simple"];

This uses the Cross-Validated 'h' to see how good the fit looks:

{fits, resids, slopes, 
          kern} = {NPFits, NPResiduals, NPSlopes, OutputKernel} /. 
          NPRegress[Transpose[data], H -> cvh, Residuals -> True, 
            Kernel -> Normal, Slopes -> True, 
            Order -> 1];

\(\(MultipleListPlot[{Transpose[data], Transpose[fits]}, 
        PlotLabel -> \*""Normal kernel, \!\(h\_cv\) = "" <> 
            ToString[cvh], 
        SymbolShape -> {PlotSymbol[Diamond, Filled -> False], 
            None}, PlotStyle -> {Hue[0.9], Thickness[ .004]}, 
        PlotJoined -> {False, True}, 
        PlotRange -> {Automatic, {\(- .4\), 3}}];\)\)

CrossValidatedH  - (con't)

Different Loss Functions

There are a wide variety of parameters that the user can fine-tune...

CrossValidatedH  - (con't)

Different Choices for LeaveOut

There are a variety of parameters that the user can fine-tune...

CrossValidatedH  - (con't)

Different Choices for Samples

There are a variety of parameters that the user can fine-tune...

CrossValidatedH - continued

Bivariate Data

CrossValidatedH - Bivariate Data - (con't)

For two- and three-dimensional data, a graph of the minimized loss function is also output so that the user can get an idea of any sorts of numerical instabilities that might be going on:

kernel = Normal; temp = Transpose[Flatten[data2, 1]]; 
    Timing[{cvh, losses, g, hrange} = {HStar, Losses, CVGraph, HRange} /. 
            CrossValidatedH[Drop[temp, -1], Last[temp], 
              LossFunction -> "MSE", Kernel -> kernel, 
              Order -> 1, LeaveOut -> 5, 
              HGridDimensions -> 16, HRange -> Automatic, 
              GraphLosses -> True, Report -> True, 
              OutputLosses -> True, Samples -> 16, 
              SearchPoints -> 4, 
              MinMethod -> "Simple"];]

CrossValidatedH - Bivariate Data - (con't)

Here is the fit that results:

MultipleKDensityPlot, SilvermanH, FastNPDensity...

MultipleKDensityPlot

?MultipleKDensityPlot

MultipleKDensityPlot quickly displays a group of distributions of a matrix of Univariate datas, color coded from first to last:

Show[
        MultipleKDensityPlot[
          N[Log[1 + Table[Age[year], {year, 1980, 2002}]]],  .4, 250, -3, 
          7], PlotRange -> {All, {0, 0.75}}, 
        AxesLabel -> {"Log[1+years]", ""}];

Bibliography

Pagan and Ullah (1999), Nonparametric Econometrics.

Silverman (1986), Density Estimation for Statistics and Data Analysis.

Härdle (1989), Applied nonparametric regression.